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Abstract 

We study the stability of the isotropic vacuum Friedmann universe in gravity theories with higher- 
order curvature terms of the form (JJaft-R"')" added to the Einstein-Hilbert Lagrangian of general relativ- 
ity on approach to an initial cosmological singularity. Earlier, we had shown that, when n = 1, a special 
isotropic vacuum solution exists which behaves like the radiation-dominated Friedmann universe and is 
stable to anisotropic and small inhomogeneous perturbations of scalar, vector and tensor type. This is 
completely different to the situation that holds in general relativity, where an isotropic initial cosmolog- 
ical singularity is unstable in vacuum and under a wide range of non-vacuum conditions. We show that 
when n 7^ 1, although a special isotropic vacuum solution found by Clifton and Barrow always exists, it 
is no longer stable when the initial singularity is approached. We find the particular stability conditions 
under the influence of tensor, vector, and scalar perturbations for general n for both solution branches. 
On approach to the initial singularity, the isotropic vacuum solution with scale factor a{t) = t^-^^ is 
found to be stable to tensor perturbations for 0.5 < n < 1.1309 and stable to vector perturbations for 
0.861425 < n <1, but is unstable as t ^ otherwise. The solution with scale factor a{t) = t^+^^ is not 
relevant to the case of an initial singularity for n > 1 and is unstable as t ^ for all n for each type of 
perturbation. 



1 Introduction 



The study of the very early universe leads us to investigate what happens to our assumptions about the 
truth of Einstein's general theory of relativity when the curvature of space and the density of matter and 
radiation approach the fundamental Planck values defined by the constants of Nature, G, c and h. The most 
natural extensions to explore as generalisations of general relativity are the higher-order theories of gravity 
that arise when the Einstein-Hilbert Lagrangian is extended by adding powers of the scalar curvature or 
the square of the Ricci tensor. As the Planck epoch is reached, or passed, on approach to a cosmological 
singularity, these higher-order terms are expected to dominate the behaviour of simple cosmological models. 
Any evaluation of what are likely initial conditions during the pre-inflationary era of a cosmological model 
should therefore be based on a full understanding of the general behaviour of cosmological models in the 
presence of higher-order gravity terms. 

Contributions to the Lagrangian from powers of the scalar curvature, i?", are conformally equivalent to 
the presence of a self-interacting scalar field and are understood [1]. In an earlier paper [2], we considered the 
effect on cosmological singularities of adding the quadratic Ricci invariant RabR""^ to the Einstein-Hilbert 
action of general relativity. The purely quadratic Lagrangian gravity theories that contain this invariant, but 
not the Einstein-Hilbert {R) term, possess an isotropic vacuum cosmological solution, in which the expansion 
scale factor, a{t), behaves as in the flat Friedmann radiation-dominated universe of general relativity, with 



a{t) = i^l"^ [3, 4]. In the case of zero spatial curvature^, this vacuum solution of the pure RabR'^'' theory 
therefore has the exact metric: 

ds^ = -df + t{dx'^ + dy^ + dz^). (1) 

Thus, we see that the higher-order Ricci stresses induce a behaviour that mimics the effect of an isotropic 
black-body radiation stress, even though no physical stress of this sort is present. Earlier studies of 
anisotropic, spatially homogeneous universes of Bianchi types I, II [4] and IX [5] showed that this special 
isotropic solution is stable against homogeneous anisotropic distortions as t ^ 0. This surprising situation is 
completely different to that encountered in general relativity (GR), when the RabR"^^ term is absent from the 
action. In GR, the expansion and 3-curvature anisotropies dominate the vacuum dynamics as t — » so as to 
produce anisotropic [6, 7], and c;ven chaotic [8], dynamics. For all pcrfcx;t fluids with pressure, p, and density, 
p, satisfying —p/3 < p < p, the isotropic solution is unstable as t and hence such isotropic solutions 
are special in GR [9]. This instability does not occur when the RabR"'^ term is present. On approach to the 
cosmological singularity, the higher-order curvature terms render the isotropic solution stable. This has all 
sorts of consequences for physical cosmology. For example, it ensures that a pre-inflationary state will likely 
be isotropic and it removes the need for the introduction of an extra physical principle, like the minimisation 
of a 'gravitational entropy' [10, 11], in order to enforce a special isotropic initial state. However, it does 
suggest that a stable state of isotropic contraction will be produced on approach to any future singularity in 
a closed universe and that may be an awkward conclusion for any theory of a gravitational entropy governed 
by its own gravitational 'Second Law'. 

The addition of quadratic Ricci terms can also create unusual evolutionary behaviour, not seen in general 
relativity. Barrow and Hervik found exact solutions which display anisotropic inflation [4, 11]. These 
solutions do not have a general-relativistic limit and are intrinsically non-linear with respect to the space- 
time curvature. 

In our first paper [2] , we extended the study of the effects of an RabR"'^ addition to the Einstein-Hilbert 
action to the situation of anisotropic and inhomogeneous cosmologies. Specifically, we investigated the 
behaviour of small scalar, vector, and tensor perturbations to the metric (1) as t 0. We found that 
there were no growing metric perturbation modes of scalar, vector, or tensor sorts as t — > 0. Thus, a small 
perturbation of the isotropic cosmological solution forms part of the general solution of the gravitational field 
equations when the RabR""^ term is present: it is an open property of the initial data space of the quadratic 
theory. 

These results immediately suggest that we should investigate whether or not the stability of isotropic 
singularities is maintained to higher order when we introduce additions to the Einstein-Hilbert action of the 
form {RabR"'^)"- We expect the situation for n 7^ 1 to be more complicated because there will no longer be 
a simple Gauss-Bonnet invariant underlying the field equations. This question of the stability of the n 7^ 1 
theories is the subject of this paper. In the absence of the Einstein-Hilbert term there is a counterpart to 
the simple isotropic vacuum solution of equation (1) in the case of general n, which was found by Clifton 
and Barrow [12]. This reduces to the solution (1) as n — > 1. It is the stability of this isotropic power-law 
solution for general n that we shall investigate. 

In section 2 we give the field equations for the gravity theory with an ii -|- A{RabR°'^)" Lagrangian and 
give the exact isotropic vacuum solutions. These solutions have two branches. We identify the physically 
interesting one that describes an expanding universe and show that as n ^ 00 the exact vacuum solution 
approaches that of a dust-filled general relativity solution with a{t) = t^/^. 

In section 3, we present the formalism for studying small tensor, vector and scalar perturbations of this 
special vacuum solution in order to determine the conditions on n for which it is stable as t ^ and the 
initial singularity is approached. In sections 4, 5, and 6 these stability analyses are carried out for tensor, 
vector, and scalar perturbation modes, respectively. The results are summarised and discussed in section 7. 
A collection of useful quantities is derived in the appendices. 

'^The Friedmann radiation solutions are also exact solutions of the pure RaiR'^^ theory in the cases of non-zero spatial 
curvature [4]. 
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2 Field Equations 

Consider a higher-order gravity theory with action 

S = j d'^x^T^ 



X 



where x, ^ and n are constants. The field equations are obtained using the general formula from Clifton and 
Barrow [12] which expresses the higher-order contributions as an additional eS^ective stress tensor: 



Gt + AP,- = ^T,% (2) 



where 



+ng^ (Y''-\cdR"' + 2(F"-i),eii'=1d + ^V'-^OR 

+{Y^-^),cR'",b + (F"-'),c5"''i«6;/ + Y^-\g'''^R,db + 2R\abR'"^)^ , 

with Y = R°-^Rab, and Gab = Rab — \Rgab is the usual Einstein tensor. 

We consider perturbations about a spatially flat, homogeneous and isotropic FRW spacetime with metric 

ds^ = -dt^ + a'^it){dx'^ + dy"^ + dz^), (3) 

with aforementioned scale factor a{t) and associated Hubble expansion rate H = -. 

In the limit where the Ricci term dominates, A ^ oo, which we expect to be appropriate in the neigh- 
bourhood of the cosmological singularity where a ^ 0, provided n > ^, the vacuum field equations reduce 
to = 0. To background order, we have: 



po = -Y''-'^{^-Y'^ + 6nY{2HH-2H'^ + 3H'^H-3H'^)+6n{n-l)Y{2HH + 3H^)j, (4) 
P° = = P^, (5) 
P« = -y"-35^ { + {"^'^ + ^^-^-^ +12H'^ + 18H'^H - 18iJ^) 

+n{n - 1)YY (sH + 36HH + 12H^^ + 2n{n - 1) (2H + 3H'^^ (^{n - 2)^"^ + YY^^ (6) 
Substituting for Y and Y in terms oi H, H, . . . gives 

po = -72(12ij"2 + SQH'^H + 36//^)"-^ jir* + QH'^H^ + IbH'^H'^ + ISH^H + 9H^ 

+n{-2HH'^H - QH^HH - 3H^H - 2H^ - 15H^H^ - 42H^H'^ - 3QH^H - 9H^) 
Wi^HH'^H + 12H^HH + 9H^H + UH'^H^ + 42H'^H'^ + 36il^ii)|. 

Hence, the Friedmann-like equation for this theory in vacuum is 

= H^H{9n^ -3n) + H^HH{12n^ -6n) + HH'^Hiin^ -2n) + H'^{l-2n) 

+H^H^{6 - 15n + 12n2) + H'^H^{lh - 42n + A2n^) + H^H{18 - 36n + 3671^) + H^{9 - 9n). (7) 
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For power-law scale factors, a = t'', and general values of n ^ l,tliis implies 



1 i P 
k = 0,k = - ± -\/3, or fc = — , 

2 6 3 

where the possible values of P are given by the two roots of a quadratic: 

3(1 - 3n + 4n2) ± ^3(-l + lOn - 5n^ - 40n^ + 48n^) 



P = P± = 



2(1 -n) 



(8) 



-4 -2 



2 in 



Figure 1: The variation of P+ with n. 




Figure 2: The variation of P_ with n. 



In the limit n ^ 1, P_ 



|, and we obtain the special a 



case studied in ref. [2]. Note also that P_ rapidly asymptotes towards 2 as n 

1 1 13 



vacuum solution of the quadratic (n = 1) 
oo, 



p_ 



2 



3n 18n2 216n3 



0(n-4), 



(9) 



and the vacuum solution rapidly approaches the behaviour of the GR dust solution with a = t^/^, see 
equation (9). P (or its real part) is greater than 3 only for the range — ^ < n < —0.390388. For the choices 
fc = 0, fc = \ ± 1^3) we must have n > 0, since 12ij^ + 36H'^H + 36il^ also vanishes. The physically 
interesting cases relevant to an initial singularity arc those with fc > 0, i.e. solutions which arc expanding to 
the future. Finally, we note that an exponential scale factor with H = constant is possible iff n = 1. 

For comparison, in a perfect fluid-filled universe with equation of state p = w/i, there is a flat FRW exact 
solution of the {RabR"'^)"' theory where the scale factor is given by 



a{t) 



(10) 



3 Inhomogeneous Perturbations 

We will now develop the formalism for studying small perturbations of the spatially flat isotropic FRW 
solutions of the {RabR°'^Y theory, which generalises the formalism developed by Noh and Hwang for the 
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quadratic (n = 1) thcjory [13, 14, 15]. We are interested in the stability of the spatially flat isotropic 
background FRW solution 



a{t)=t^, (11) 

where 

_ 3(1 - 3n + 4n2) ± V3(-l + lOn - 5n^ - 40n^ + 48n4) 

The general perturbed metric may be written as 

ds^ = -a'^{l + 2a)dr]^ -a^Badrjdx" +0^(50,(3 + Ca0)dx"dx'^, (13) 

where is a conformal time coordinate that is related to the comoving proper time, t, by dt = adrj. We can 
decompose the perturbation variables into their scalar, vector and tensor parts in the standard way, as in 
[2], by writing 

S„ = 2/3„ + 2S„, 

Cal3 = 'i(t>5al3 + 27 + 2C(„,/3) + 2Ca/3. 

There are four scalar perturbation variables, a.p.cp and 7, two vector variables, Ba and C'a, and one 
tensor, Caf3- The quantities and Ca arc divergence-free, i.e. B"^a = = C"q,, and Caf3 is transverse 
and trace-free. These three typos of perturbation evolve independently of each other at linear order. We 
will determine the equations which describe their time evolution and then solve each of them to determine 
whether the metric perturbations to the special solution arc stable as t ^ i). In the n = I case the problem, 
the equations, and their solutions will reduce to those of [2]. In this way wo establish the ranges of n values 
for which the special isotropic vacuum solution is a stable initial condition for the higher-order theory. 

4 Tensor (gravitational-wave) perturbations 

The expansion of the metric around the spatially flat Friedmann solution now takes the form 

ds^ = -dt^ + a^{6a0 + 2Ca0)dx°'dx^. 
The tensor Ca^ is trace-free and transverse, i.e. 

CS = = C^,a (14) 

and C = C(x, t) . 

The n = 1 case was solved exactly in [2] for perturbations about a{t) = t^/'^. Here, we want to perturb 
an isotropic background solution which has a{t) = t^^^ with 



_ 3(1 - 3n -h 4n2) ± v'3(48n4 - 40n^ - Sn^ + lOn - 1) 

The important quantities to linear order in the perturbation are given in Appendix A. In the limit where 
the higher-order terms dominate, the perturbed fleld equation is 

6P^ = -nY'^-'^(^C^+6HC0+SH'^C^-{3H + 21HH+18H^)C^-2^C^-2H^C^ 

+ (4iJ + 8H^ + ^) - n{Y-'y (^2C$ + 9HC^ + 3H^C^ - 2^C^ + H^C^^ 

-n{Y--'Y (c^ + 3HC^ - ^C^) . (16) 
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4.1 Lcirge Scales 

In the long-wavclcngth limit, on super-horizon scales, we can neglect terms involving AC, A^C and AC. 
For a{t) = t^/^, we have H = ^ and Y{t) oc so the equation for the perturbations becomes 



= SP^ = -nV'-U C$ + {-8n + 8 + 2P)-^+ 1 16n'^ -28n+ 12 + UP- 12Pn+—j 



and so 



C oc 

= A^A^ + (-8n-t-2 + 2P)A^ + ^-4n- 1 + 16n2 + 6P(1 - 2n) + A 

-2 + 12n - 16n2 + P(2 - 16n + IGn^) + p2 _ _ 2p3^ (^g) 
The four roots of this are 

A = 0, Ai, A±, (19) 



where 



A± = i(Ai±v^), (20) 
Ai = -l-P + 4n (21) 
A2 = ^P2 - 14P + 8Pn + 9 - 24n + IGn^ (22) 



464n4 - 392n3 + 9n2 + 64n - 13 ± (36n2 - 1 In - 3) v'3(48n4 - 40n^ - 5n^ + lOn - 1) 

2(1 -n)2 ■ 

Ai and A2 are real whenever P is real, i.e. n ^ (—0.47942,0.110873) 2. For A± to be real we need A2 > 0. 

We are interested in the signs of the possible values of A in order to determine the behaviour of 
gravitational-wave perturbations of the isotropic solution as t ^ 0. If any 3?(Ai) < 0, the solution is 
unstable as t ^ 0. Otherwise, we need to look at the stability problem to second order, due to the presence 
of the zero eigenvalue. For reference, we recall that for the n = 1 theory, studied earlier [2] , we had 

C(x, t)(xa + +^t + 5t^''^ 

and there were no diverging metric perturbation modes as t ^ 0. Let us now analyse the situation in the 
more complicated n 7^ 1 case: 

4.2 Solutions with P = P+ 

First consider the case P = P+, so that 



3(1 - 3n + 4n2) + ^3(-l + lOn - 5n2 - 40n3 + 4871^) 
~ 2(1 - n) ' 



^The values here and in the tables that follow are approximate numerical roots of the appropriate polynomials. 
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the stability is decided by the quantities 

-20n2 + 19n - 5 - + lOn - 5n2 - AOn^ + 48n^) 

2(1 - n) ' 
46An^ - 392n3 + Qn^ + 64n - 13 + {36n^ - lln - 3)y3(-l + lOn - Sn^ - 40n3 + 48n4) 

2(1 -n)2 

u 

. I . . . n. 



Figure 3: Power-law exponents, Ai, A±, versus n for tensor perturbations with P = 





Values 
P+ 


of Xi for diffe 
Ai 


jrent v 

A2 


alucs of n 7^ 1, taking P = 
A+ A_ 


P+ 

Remarks 


n < 
-1.30084 


M 

> 


M 

< 


K 
> 


R 

> 


R 
< 


Unstable as t — > 0. 


-1.30084 

<n< 
-0.47942 


M 

> 


R 
<0 


K 
> 


R 

< 


R 
<0 


Unstable as f — > 0. 


-0.47942 

< n < 
0.110873 


c 

3fi(P+) > 


c 

3?(Ai) < 


c 


c 


c 


Unstable as t — > 0. 


0.110873 

<n< 
0.452692 


M 

> 


M 
<0 


K 
< 


c 

5R(A+) < 


c 

^{X-) < 


Unstable as t — > 0. 


0.452692 

<n< 
0.5 


M 

> 


M 
<0 


K 
> 


M 

< 


M 
<0 


Unstable as f — > 0. 


0.5 
<n< 
1 


M 

> 


M 
<0 


R 
> 


M 

> 


M 
<0 


Unstable as f — > 0. 


n > 1 


M 

< 


M 
>0 


M 

> 


M 

> 


R 
<0 


Universe contracts. 

P+ < 



We note that for the solutions with a{t) = t^+/^, there is always a negative eigenvalue, so are unstable 
for any n as t — > 0. For n > 1, P+ < 0, so this corresponds to a contracting universe, in which we are not 
interested here. However, we have 3f?(P+) > 3 for 1 > n > ^ and n < so we need to be careful that 

the instability for these n is not arising from the negative curvature contribution characteristic of the Milne 
universe. We expect the overall assumption that the higher-order Ricci terms dominate the GR terms in the 
neighbourhood of the initial cosmological singularity to hold so long as n > 1/2. 
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4.3 Solutions with P = P_ 



Now consider the second case, with P = P_, which turns out to be the most physically relevant for consid- 
eration of the effects of higher-order (n > 1) corrections. Wc have 

3(1 - 3n -I- 4n2) - v'3(-l + lOn - 5n2 - AOn^ + 48n4) 
~ 2(1 -n) ' 

with the stability decided by 



Ai = 
A2 = 



-20n^ -I- 19n - 5 -f a/3(-1 -|- lOn - Sn^ - 40n^ + 48n4) 
2(1 - n) ' 
464n* - 392n3 -|- Qn^ + 64n - 13 - {36n^ - lln - 3) V3(-l + lOn - Sn^ - 40n3 48n4) 

2(1 -n)2 ■ 




Figure 4: Power- law exponents, Ai, A±, versus n for tensor perturbations with P = P- 
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Values c 
P_ 


)f Xi for difFe 
Ai 


rent vt 
A2 


dues of n ^ 1 
A+ 


, taking P = 
A_ 


Remarks 


n < 
-0.47942 


K 

> 


R 

< 


R 

> 


R 

< 


R 

< 


Unstable as t — > 0. 


-0.47942 

<n< 
0.110873 


c 

3?(P_) > 


c 

3?(Ai) < 


c 


c 


c 


Unstable as t — > 0. 


0.110873 

<n < 
0.159452 


R 

> 


R 

< 


M 

< 


c 

3fJ(A+) < 


c 

3?(A_) < 


Unstable as t — > 0. 


0.159452 

< n < 
0.169938 


K 
> 


R 

< 


M 

> 


M 

< 


R 

< 


Unstable as t 0. 


0.169938 
<n < 
0.25 


K 

> 


R 

< 


M 

> 


M 

> 


R 
< 


Unstable as t ^ 0. 


0.25 
<n < 
0.5 


R 

< 


R 

> 


M 

> 


M 

> 


R 

< 


Universe contracts. 

P_ < 


0.5 
<n < 
0.520752 


R 

> 


R 

> 


M 

> 


K 

> 


R 

> 


Stable as i ^ 0. 


0.520752 

<n < 
0.989666 


R 

> 


R 

> 


M 

< 


c 

3?(A+) > 


c 

3?(A_) > 


Stable as t ^ 0. 


0.989666 
<n < 
1.1309 


R 

> 


R 

> 


M 

> 


K 

> 


R 

> 


Stable as t ^ 0. 


n > 1.1309 


R 

> 


R 

> 


K 

> 


M 

> 


R 

< 


Unstable as t ^ 0. 



Wc saw that for n > 1, P+ < 0, so for an expanding universe, the only relevant value for P with n > 1 
is P = P_ . The first table shows that this can only be stable as t ^ if 

^ < n < V2 cos 

In particular, for all integers n > 1, the exact isotropic solution with a{t) = t^-/^ is not a past attractor 
as t ^ 0. Thus, it appears that the quadratic (n = 1) case studied earlier was exceptional and the stability of 
the isotropic singularity found for that case does not extend to higher-order corrections to general relativity 
with n > 1. 

5 Vector (vortical) perturbations 

We have shown that for gravitational-wave perturbations there is a very small range of values of n, given 
in eqn. (23), for which the perturbations are stable as t — > 0. We now consider the vortical perturbations, 
which are of vector-type. The metric is: 

ds^ = -dt^ - 2aBadtdx°' + a^iSaff + 2C(^a,f3))dx°'dx'^ , 



- arccos 
3 



-V2' 



1.1309. 



(23) 
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where = = C"„. 

The energy- momentum tensor is decomposed as usual [16], 

Tab = IJ.UaUb + phab + QaUb + QbUa + T^ab, (24) 

where hab = Oab + UaUb, QaU^ = = TTabu'' and TT^ = 0. The fluid four-velocity Ua and the energy flux qa are 
decomposed as (using t as the time variable, index '0' denotes t): 

?i" = 1, uo = -1, 

qo = 0, qa = aQa- (25) 

The energy-momentum tensor is then 

= (m + P) (Vo. + - B^) , ST^ = n^. (26) 

We can also decompose the perturbation variables as 

B„(x,i) = 6(t)F„(x), a = cy„, Ay„ = -fc2y„, 

Vo, = vY^, Qa = qYa, ni=pnTY^ (27) 
and introduce the gauge-invariant variables [14] : 

V„ + -^-B^=(v+^—-h\Yo, = v^Y^, (28) 



IJ' + P V M + P 

K. + -^+C^,= (z; + ^ + c')y„ = v^Y„, (29) 



So + C = (W. - l^^)ya = *ya, (30) 

where a prime denotes a derivative with respect to the conformal time variable t]; Vu and Va may be 
interpreted as the velocity variables related to the vorticity and the shear respectively. 

Wc will work in the "C-gaugc", i.e. we act Ca = 0, which completely fixes the gauge condition. Then, 
using the quantities presented in Appendix B, we find that the perturbed parts of the tensor Pj^ are: 

6P^ = 0, (31) 
+n{Y^-^y I^B^ + H-bA, (32) 



\ 2a , 

6P^ = -ny^-iijB^",^) +3ifS(",^) - (^3H + 6H^ + ^^B^''j,)-A{H + 6HH + 4H^)B^''^p)^ 

-n(y"-i) -i +2ifi3(",^)}. (33) 
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Combining the equations (26)- (33), we have 



2a 



+ A 



A 
2a' 



2 



A2_ 
2a3 



(34) 
(35) 



and so for vanishing anisotropic pressure of the matter part, pwT = 0, angular momentum is conserved 
exactly as in the quadratic case and a^T^ = flYa{x) is a constant in time. 



For a = t^/^, if = ^, and y = 1^ (3 - 3P + P^) oc r", so we have therefore: 



po = ny"-l- 



2^3- 



4 + 



4n B^ 



f2 



4P- 



8P2 8Pn\ 



2P 

t 3 



(36) 



5.1 The oo limit 

In the limit A ^ oo, where the GR term can be neglected and the higher-order Ricci terms dominate, we 
have 



4P2(3-3P-(-P2) 



1 



| + 7!l(-4(n-l) + -)* 



+ 



^2 



2 2t 

APn 6P 4P2 



2t¥ 



(37) 



Defining the constant 



9. 



2Q 

ni^ V4f^(3-3P + P2) 



we have 



9.t 



An— A- 



1 



* + J ( -4(n - 1) + - ) + _ ( -4Pn + 6P - — ) * + 



4P^ 



If we take the long-wavelength limit, i.e. we drop the last term on the right-hand side, then we have to solve 



m4»-2-¥ 



P 



= r^f + f _4(n - 1) + ^ * + o -'i^'" + 5- 



4P^ 



The complementary function (LHS = 0) is solved \sy = t^, where 



= e 



3-4n + 



P 



-4Pn + 6P 



4P- 



^ = ^± = ^ (-9 - P -M2n ± V3(27 - 72n 48n2 - 42P + 24nP + 11P2)^ , 

= ^(Ai+8(n-l))±iv^, 

where the were defined in (21) and (22), and for stability as t — » 0, we need 3?(C±) > 0. The additional 
mode from the particular solution has ^ ~ /;4n-2-2£ _ ^|(Ai+2(n-i))_ signs of these exponents for 

different values of n are summarised in the tables which follow and their values are plotted in figures 5 and 
6 for P+ and P_ respectively. 
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5.1.1 Solutions with P = P+ 



For P = P^, wc have 

> for n > i and n < -if, whilst 3?(C_) < for ah n. Finally, from the 
particular solution, the exponent 4n — 2 — ^ is positive for n > 1 and negative (and hence unstable as 
f — > 0) for n < 1. Thus, for the solution branch P = P+, the vector modes are unstable as t — > for all 
values of n. 



40 




Jl, 


-i -1 

--30 

^^y^ -40 




1 4 



Figure 5: Power-law exponents, 4n — 2 — ^_|_, versus n for vector perturbations with P = P+. 





Value 


s of for dil 

u 


Ferent values 

i- 


of n ^ 1, taking P = 1 
4n-2-^ 


3^ 

Remarks 


n < — ^ 


M 

> 


K 

< 


M 

< 


R 

< 


Unstable as t ^ 0. 


Vi 
18 

< n < 
-0.47942 


K 
> 


R 
> 


R 
< 


R 

< 


Unstable as t ^ 0. 


-0.47942 

< n < 
0.110873 


c 

5R(P+) > 


c 

< 


c 

< 


c 

5?(4n-2 - ^) <0 


Unstable as t ^ 0. 


0.110873 

< ri < 
0.5 


R 
> 


R 

< 


R 
<0 


R 

< 


Unstable as t ^ 0. 


0.5 

< n < 
1 


R 
> 


R 

> 


R 

< 


R 

< 


Unstable as t — > 0. 


1 < n 


R 
< 


R 
> 


R 
<0 


R 

> 


Universe contracts. 

P+ < 



5.1.2 Solutions with P = P_ 



For P = P-, 5R(C+) > for i < n < i and n > 0.861425, 3R(C_) > for 0.861425 < n < 1, whilst 4n-2 - ^ 
is positive for n > 0.5 and negative for n < 0.5. For this branch, the vector perturbations are stable to linear 
order as t ^ for 



1 > n> 



1 

36 



25-|-2V23sinh 



larcsinhf^ 
3 V23\^ 



0.861425 



(38) 



and are unstable for all other n. 
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Figure 6: Power-law exponents, 4n — 2 — ^_|_, versus n for vector perturbations with P = P- 
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4n - 2 - — 


Remarks 


n < 
-0.47942 


M 

> 


M 

< 


R 
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R 

< 


Unstable as t — > 0. 


-0.47942 

<n< 
0.110873 


c 

3?(P_) > 


c 

m+) < 


c 

m^) < 


c 

3?(4n-2-2£) <0 


Unstable as t — > 0. 


0.110873 

< n < 
0.159452 


M 

> 


c 

5R(?+) < 


c 

3^(4-) < 


M 

< 


Unstable as t — > 0. 


0.159452 
< n < 
0.25 


M 

> 


R 

< 


R 

< 


R 

< 


Unstable as t ^ 0. 


0.25 

< n < 
0.5 


M 

< 


R 

> 


R 

< 


R 

< 


Universe contracts. 

P- < 


0.5 
< n < 
0.520752 


M 

> 


R 
< 


R 
< 


R 

> 


Unstable as t ^ 0. 


0.520752 

<n< 
0.861425 


M 

> 


c 

m+) < 


c 


R 

> 


Unstable as i ^ 0. 


0.861425 

< n < 
0.989666 


M 

> 


c 

m+) > 


c 

3?(e_) > 


R 

> 


Stable as t ^ 0. 


0.989666 
< n < 
1 


R 

> 


R 

> 


R 

> 


R 

> 


Stable as i ^ 0. 


1 < n 


R 

> 


R 

> 


R 

< 


R 

> 


Unstable as t ^ 0. 
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6 Scalar perturbations 

We will now consider scalar perturbations. The metric for the general scalar-type perturbation takes the 
form 

ds^ = -{l + 2a)dt^ -2a(3^^dtdx°' +0? {6ocf){l + 2(i)) + 2'y ,^lj)dx°'dx'^. (39) 

We use the proper time, t, as the time variable and also define the quantities, x = + 07), / = y"~^ = 
{RabR""^)^'^ ■ We use overbars and deltas to denote background and perturbed quantities, so that in general, 
A = A + 6 A, and in particular / = ^Sf = {n — SY. The important quantities to linear order in 

the perturbation are given in Appendix C. 

Using the gravitational field equations, we obtain the complete set of gauge-ready equations for the per- 
turbed variables: 



Energy: 



nf 



--6Y + 12H^a + {12HH + 18H^)a + 12{AHH - AH"^ + 6H^H - m^)a 



-{IQH + MH'^)^a - 2— a - 12H^ + {2AH - im^U - (12H + 36HH - 72H^U 

, A A2 A2 ■ 

-{4H + 20HH - 20H^)^x + 2^X + '^^^x 



+nf 



UH^a + 2AH{2H + 3i?2)a + ^H^a - 12if^ - {12H + 54if2)(/) + ^H^(, 



-AH^X-{^H + WH^)^X 



+n 

Momentum: 
Pa = nj 



{-12HH + 12H^ - ISH-'H + 18H'')df - 3if (4if + 6H'^)6f + {AH + m'^)^5f 



t2 tj 



A 



-4Ha - (12H + 12H^)a + (-12H - 4SHH)a - 2^a + 4H^a 



+40 + \2B(^ + (24if)0 - 2-2 ( 



%H%^ + 2%x - 2i?4x + (2if - 8i?2)4x 



+nf 



-AHa - {12H + im^)a - 2^a + 40 + 12i/0 - 2^(l> + 2^x + 2i?4x 
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Trace: 



nf 
+nf 



A 



A 



A 



.A 



12Ha + 3{16H + 2AH^)a + 4— a - 120 - 36if(/) + 4^(/) - 4— x - '^H—x 



A 



.A 



2AHa + 3(28i/ + 42if^)d + 3(32ff + lUHH + 48if'^)a + 8— d - 127J— a 



-2A4> - 108H'<j) - 108(H + H^U + 16^^ - 16H^(t> 



i^X - ^H^x - (20ir - AH^)^x 



12i?d + 3(16ij + 2AH'^)a + 3(24ff + IQSHH + 36i?3)d 



A A A A^ 

+3(16ir + Q%HH + 72iJ^ + lUH'^ma + 4^d - 16il^d - (20H + 34ij2) « - 2— 

-12'^" - 72H'^ - 3(36i/ + 36H^y(j) - 3(24ff + 72HH)ij) 

A- A . , . o.A A2 

+16^./) + WH^^ - {16H + 8H^)^(j) - Q^(t> 



-44x - (12ir - l2H^)^x - i^QH - 8H^)-x + 2—x + 6H—x 



A 



A A2 
X + 2— 



a4 



-3n 



Sf{4H + 6if2) + ^/(8# + 36HH + 12H^) + df{4H + 2AHH + 12ij"2 + im'^H - 18H'^ 



+4n{H + 3H^)^Sf 



A 



Trace- free: 



Pia tra pa 

13 ~ 3"/3-^a 



^ (v-V^ - ^<5^AJ f[a + ct>-x-Hx\ 
+7 [2d + 3ifa -2(1)- H(l)-2x- Hx+ {AH + 7H^)x 

+7 



a + 5Ha + {AH + 8H^)a + - 50 - 17if0 - {AH + 8H^)(f) + 3^( 



-X + {6H + 9H^)x + {5H + 21HH + 8H^)x -^X- 3H^x 



-2{H + 3H'^)5f 



with 



SY 



-12H{2H + 3H'^)a - A8{H'^ + m'^H + 3H'^)a - A{2H + 3H'^)^a + 12{2H + 2,H'^)'4) 



+72H{H + 2H'^)4> - 8{H + 3if2) ^ + /^{2H + 3H^)^x + &H{H + 3H^)^x 



(40) 



The special unperturbed solution has a = = £ and F = (3 - 3P + P^) cx t''^. 
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6.1 Lineeirising about the special solution in the zero-sheeir gauge 

We now take the large-scale limit and c;hoose the zero-shear gauge (x = 0) and linearise about the special 
flat FRW solution with a = t^l^ . For P = P±, the equations simplify to: 
Energy: 



'Ap2 



8P'' 



(-6 + 6P - + n(12 - VIP + 3P^)) \ x 



jia + (2 - 4n + P)a; - -| {t^"<i) + (5 - 4n + P)t^ + (4 - 8n + 1P)<^ | 



Momentum: 
<5P° 



/4P2 

n(9^(3-3P + P^) 



8P3 
27*6 



(6 - 6P + p2 - n(12 - 12P + 3P^)) 



V„ jta + (1 - 4n + P)q! - {t^"<i) + (4 - 4n + P)4 + (2 - 8n + 2P)<?i) | 



Trace: 



4P2 



n-2 



8P3 

9*7 



(-6 -h 6P - P^ + 3n(P - 2)2) 



t^a + (4 - 8n + 2P)fa. + (2 - 4n + P)(l - 4n + P)a. 



3^ 
P 



(t^'^' + (8 - 8n + 2P)t^"^ + (7 - 4n + P)(2 - 4n + P)i^ + 2(2 - 4n + P)(l - 4n + P)^ 



Trace-free propagation: 



1 /4P2 \"" / 1 \ 4P2 

n^(^(3-3P + p2)) (v"V,--.,"a) |^(3-3P + P2)a 

-(72 - 69P + 27P2 - P^ + 6n(-12 + 14P - IP"^ + P^))a 
-(3 - 3P -f p2)(27 -h 12P - 4P2 -t- 3n(-21 - 5P 2p2) -h mn^)a 



4p2 

'27*5' 
16P2 

'"81*6" 
4P2 



+^(-3 - 3P + P2 + 6n(-2 + 3P - P'^))<^ 



4p2 

'27*5 
16P2 



(-72 ^ 57P - 33P2 7P3 12n(2 - P)(3 - 3P 2p2))</) 



+^^^(3 - 3P + p2)(27 - 2P2 + n(-63 + 3P) + m-a^)^ 



The energy and momentum equations together imply that 

3 

a = —(*?!' + ?:> + ao), 



(41) 



where ao is a free constant. In addition, this also satisfies the trace equation. 
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Finally, we use the trace-free equation to calculate that 



a 
Pi 



= I [(1 + Pi)<t>it'' + (1 + P+)^+t'+ + (1 + P-)^-t'- 
27 - 2P2 + „(_63 + 3P) + 36n2 
~ 27 + 12P - 4P2 + 3n(-21 - 5P + 2P2) + 36n2 
P 



= -3 + 4n - 



3 

= i (-15 + 12n + P ± ^3(27 - 42P + 11P2 - 72n + 24nP + A%v? fj 



^± + IT - 1. 



(42) 

(43) 
(44) 

(45) 



where (/'oi'/'i and (/)_|_ are free constants, and ^_|_ were defined in (38). We note that these are the same 
power-law exponents as for the vorticity perturbations plus -y — 1. These exponents are shown in figures 7 
and 8. For all n, and either choice of P = P±, at least one of these exponents has negative real part, and 
hence the isotropic vacuum solution is unstable as f — > 0. 




Figure 7: Power-law exponents, pi,p^, versus n for scalar perturbations with P = P_|_. 
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Figure 8: Power-law exponents, Pi,p±, versus n for scalar perturbations with P = 



7 Summary 

We can now summarise the results for the linearised tensor, vector and scalar perturbations about the 
spatially flat vacuum FRW solution (11) of the theory with Lagrangian {RabR°^^Y ■ The general perturbed 
metric in the neighbourhood of the isotropic vacuum solution a — t^^^, given by equation (11), has the form 

rfs2 = _(i _|_ 2a)dt^ - aBadtdx" + a^{SafS + Cap)dx°'dx^, (46) 
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and the perturbation variables may be decomposed into their scalar, vector and tensor parts by writing 

Ba = 2/3 + 2Ba, 

Cai3 = '^<t>^al3 + "^1 ,al3 + 2C(a,/3) + 2(7^/3. 

In the gauge defined by /? = = 7 and Cq, = 0, the general solution of the linearised equations is given by 

(47) 

(48) 

(49) 

-i+f , (50) 



= a^(x)+i^i6^(x)+t^+c^(x)+t^-rf^(x) 



Br, = 



t«+y«(x) + t€-yi^)(x)- 



a = -aocpo - 



0i(x)r3+^"-^ 

6- 12n + P 



2(3 - 6P + P2 ^ 6n + 6nP) 
^+(x)t«+-i+* +(/)_(x)i«--i+-; 



-3+4n- 



^ + ( % + 1 ) <^+*«--^^- + 



3^- 



+ 1 



where 



Ai 
A2 

A± 
^± 



-l-P + 4n 



14P+Hp2 



24n + 8Pn + 16n^ 



2v.^±v/A. 
i(-9-P-12n)±i 
^(Ai+8(n-l))± J 



27 - 2P2 + „(_63 + 3P) + 36n2 



36n2 



27 + 12P - 4P2 + 3n(-21 - 5P + 2P2) 

For the solution branch defined by P = P+, each type of perturbation (tensor, vector and scalar) is 
unstable as t — > for all values of n. For P = which is the only physically relevant value of P for n > 1, 
the tensor perturbations are stable to linear order as t ^ for 



1 



< n 



< V2 



cos 



1 



arccos 



and the vector perturbations are stable to linear order as t 



1 > n > 



1 

36 



25 + 2\/23sinh 



-arcsinh 




1.1309 



0.861425. 



,23V23, 

0. The scalar perturbations are unstable as t ^ 



For all other n, these perturbations are unstable as t 
for all n. 

In conclusion, in our earlier work [2] we discovered that isotropic cosmological models in theories of 
gravity formed with a quadratic Ricci term added to the Einstein-Hilbert action are stable on approach to 
an initial 'Big Bang' singularity. This is quite different to the behaviour of general relativistic cosmological 
models, where isotropy is strongly unstable in this limit in vacuum [6, 7]. In this paper, we have analysed 
the more complicated problem of cosmological evolution in the presence of arbitrary powers of the Ricci 
term in the Lagrangian. We have found that the behaviour displayed in the quadratic case was special. 
Isotropic power-law solutions of the sort found by Barrow and Clifton [12, 17] still exist in vacuum and with 
a perfect fluid when a term proportional to (Pa^P"'')" is added to the Einstein-Hilbert action for general 
n ^ 1. However, both solution branches of these special isotropic solutions are unstable to the growth of 
small metric perturbations as i ^ 0, and so the quadratic case with n = 1, in which these perturbations are 
bounded in this limit, is special. 

Acknowledgements: Jonathan Middlcton acknowledges a PPARC/STFC studentship. We would like 
to thank Timothy Clifton, Sigbj0rn Hervik and Kei-ichi Maeda for helpful discussions. 
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A Relevant tensor quantities using t as the time variable 

Using t as the time variable and H = ^, the important quantities to linear order in the gravitational 
wave-type perturbation are: 



i a/3 

rQ 
0/9 

ra 

r 

g I- cd 

Rq 



— pa Tfb 



R'"^R^ 
R'"^R 



cd/3 



URl 
URl 



UR 



a''[H{5o,p + 2Co,f))+Co,p] 

HSfj + Cf, 

^ /3,7 + 7,/3 — 0^/3, 

otherwise 

3(ij"-hi?^) 
i?° =0 



ttfi — 



R = 



cdO — 



6{H + 2H^) 

12(ij"2 + 3HH^ + 3H^) 

-3{H + H'^){H + 3H'^) 

= R'"'R%dQ 

-(3ij2 + 8HH^ + 9ff*)(5^ - {3H + 2H^)C^ - {7H + 6H^)HC^ 



= a 



3{H + 2HH 
= Voi?° = 

{ii + 6HH)S^ + C$+ 3HC^ + 3HC1 - + 2H^C^ 

A 
7a 



-2HHSa0 + HCa0 + (Bil^ - 2H)Ca(i -H[4H 



Cap 



6HH 

3{H + AHH) 


-3(^ -I- bHH + 2H'^ + 2H'^H) 
= URZ 

-{H + 9HH + 6H^ + 22H^H)5"p 
(7| + 3HC'^ - ^C1 



D1 + 3HD1 -[2H' + ^ 



-6 {h + 7HH + 4ij^ + 12H^H^ 
and for any scalar function f{t) of time only, it holds that 



fit);°o = -f 
Of it) = -f-3Hf 



8HHC% 
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Thus, since Y = Y{t), 

□ = {l-n)Y''-^ (^{n-2)Y'^ + YY + 3HYYy 

B Relevant vector quantities using t as the time variable 

Using t as the time variable, -ff = f and -B'",/3) = | + -B^'"), the important vector quantities to Unear 

order in the perturbation are: 



rgo = 

= —aHBa 

= 

i?o = 3{H + H^) 



2a 



i? = 6{H + 2H'^) 
Y = RIR^ ^ 12{H^ + 'SHH^ + 'SH^) 
R^'^R^M = -3{H + H^){H + 3H^) 

r>cd pO tt2 a p 

a 
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Vc-Rq 



3(if + 2//^) 

2a a 



= 0{B) 



-2a'^HH5ci3 + a 



HB, 



(«,/3) 



2{H - H'^)B(^a,0) - 



a7 



A ^ 



VpRi = —B^-aH[2H + 



2a 




□ii^ = -3(^ + 5HH + 2H^ + 2H'^H) 



t2tt 



A .. A . . o A 

—Ba + H—Ba - 2{H + H^)-B, 
2a 2a a 

-{H + 9HH + 6H^ + 22H'^H)5'i - 



rBa 



2a3' 



+ 



and for any scalar function f[t) of time only, 

f{t);°0 = -f 

f{t).\ = 

Note that i? = R{t) and y = Y{t) are the same as in the gravitational wave case. 
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C Relevant scalar perturbation quantities 

Taking t as the time variable, the metric takes the form 

rfs^ = -(1 + 2a)dt'^ - 2aP^^dtdx°' + a^it) {S^fjil + 2(1)) + 27 dx°'dx^ 
The important scalar quantities to linear order in the perturbation (using H = ^, x = a{(3 + a-y)) are: 



^ 00 

^ a/3 



■ 00 



g^'^ro^ = 3H-a-6Ha + 3(j) + 



Rq 



= a 



{a - aHI3)^, 



5ccp{H-2Ha + 2H<j) + i^^ + {2H^+^^ ^ 



{a-aHf3-af3)^' 



pa 

^0/3 — 



Oc 



3i? + d + 3(/) + A7 
(a + 3.^ + A7),„ 



A 



cd 



1 



(-a + a/3 + 2aif/3 - + A7) 



|a 



= 3(ij" + _ 3ijQ, _ 6(ij- + H'^)a - 4a + 3^ + 6il(?i + 4x 
= -2{Ha-'(P),oc 

a'' 

= (^H + 3H^ -Ha- 2{H + 3H^)a + ^ + QH^) - ^f/* + ^^xj 



R = 6{H + 2H^)-6Ha-12{H + 2H'^)a-2—a + 6'<j) + 24H^-4- 

A A 

+2—x + 4H^X 



Y 

R'^'^RPccLO 



rtcd pO 
^ ^ cda 



12(ij"2 + 3H'^H + 3il*) - 12H{2H + 3H'^)a - A8{H'^ + 3H'^H + 3H'^)a 
-A{2H + 3H'^)^a + 12(2ij" + 3H'^)^ + 12H{H + 2il2)^ _ 8(ij- + 3H'^) ^ 



A 



A 



-A{2H + 3H^)—x + SH{H + 3H^)—x 



A 



F + (5r with 

12(ij2 + 3i?2^ + 3i?*) 

-3(7?^ + 4iJ^ij" + 3ff^) + 6ff (ii + 2if2)a + \2{ii^ + 4il2ij + ZH'^)a 
+2{H + 2i?^)^a - 6(i? + 2il2)^ _ 12H{2H + 3H'^)^ + 4{H + H'^)^(i 

-2{H + 2H')^x - ^H{H + H')^x 
AH'^{Ha - (/)),„ 
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Vo-Ra 

Vo-Rq 



A 



A 



-{H + 2H^)^x - {HH + m^)^x 



1 



(3^ + 2H'^)a + {H + 2H'^)(j) - {3H + 2H'^)x - H{H + 2H'^)x 



|a 



3(F + 2HH) - 3Ha - (9H + 6H^)a - (i(H + 2HH)a - + 2H^a 

+3'</)' + 6H',j) + 6H(P + 4x - 2H^x 
-2Ha + (2ij2 - 4H)a + 2^ - 2H4) + 2aHH0 



2Ha + {AH - 2H^)a - 2^ + 2H^ + 2a{H - HH)P 



I" 



{H + 6HH)6'^ + S'^ 



-Ha - 3{H + 2H'^)a - 2{H + 6HH)a 



+■(/.■ + 6H'^ + QH^ -^^ + ^H^cj) + H^x +{H- 2H^)^x 



1 



-a + 2Ha - ^ + 2H<I> + x - Hx+ {H - 2H^)x 



Va-Ro 



Vc-Rq 



A •• - A 

-SHa - {6H + 2H^)a - + 3(?i + 2H(j} +^x + 2aHH(3 



aH. 



a/3 



-2HH + 2H^a + 8HHa + H-^a - 2H(p - 2H(j) - 4HH(f> - H-^<. 



A .A ■ 
-H^X + H^^X 



-3Ha + 2'(t>- H(t> + Hx + {H^ - 2H)x - Aa^HH-f 



,al3 



2HH - 2H'^a - AH Ha - H^a + 2H'^ + 2£r^ + H^(j> + H^x - H'^^X 



3Ha -2(j) + H(l)- Hx + {2H - H^)x 

2 



1/3 



3{H + AHH) - 3Ha - {9H + 12H'^)a - 6{H + AHH)a - + ^-H"^" 
+3^ + 12H'4> + 12H^ - 2^^ + AH^cf) + + (2^^ - 4Jf')4x 

-Ha - 2{H + 3H^)a + + 6i?^ - + H^x 
-2H5p^{Ha - </. + oi?/?)!" - 2H5'^{Ha - 4>),0 + ^{-a - + x + ^^x)'" 
-3Ha - 6{H + 2H^)a -^a + 3^+ 12H^ - 2^(j) + ^X + ^H^x 



Pi 
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3{H + 5HH + 2H^ + 2H'^H) + 3Ha + {12H + 15H'^)a + {18H + 57HH + 6H''^)a 



t2 n 



+12{H + 5HH + 2H'^ + 2H'^H)a 



-^a - AH^a - (8H + AH'')^a -ra - 3 , 

a'' a'' a* 



A . o A .. . , A A^ 

+H^X + m + 3H^)^x - i^H - 2HH + 8H^)-x + — 

2Ha + 6{H + H^)a + {6H + 12HH - 16H^)a - 2H^a 



-2j - 6H'^ + 16H^^ + 2^4) - + {AH^ - 2H)^x 



X 



= -(h + 9HH + 6H^ + 22H'^H)6'^ + 6] 



r2 ir\xa 



Ha + {AH + 9iJ2)a + (6H + 39HH + 22H^)a 



+(4ij + 3&HH + 24ij"2 + 88H'^H)a - H^a - {2H + AH^)^a - - QHJ 



-{12H + 22H^U - (9H + AAHHU + 2^'^ + hH^i() - {2H + AH"^)^, 

a'' a'' 

-H%x - {^ii + ^f')4x - (2if + 3HH - AH^)^x + H^X 



A2 



+— 



a-Ha-{2H+ 12H^)a - + ^ + 7H^ - {2H + AH'^)(j> - 



-X + 5H^X -{H + 5HH - AH^)x +^X + H^x 



OR = -6H-A2HH-2AH'^-72H^H + 6Ha + {2AH + A2H'^)a + {36H+17AHH+72H^)a 

... . A A ■ A A^ 

+24(ir + 7HH + AH^ + l2H^H)a + 2— a - 8H—a - {16H + 28H^) — a - 2— a - 6 (j) 

(X a cii (X 

-A2H4> - 2A(2H + 3H^)6 - (A2H + lAAHHU + 10^6 + 2QH^6 - 8(H + H^)^^ - 4^< 

a-' a* 

-2^X - 2H^x - (4ij - 8H^)^x " i^OH + 12HH - 8H')^x + 2^X + 4ff^x 
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and, for a general scalar function f{x,t), 



/;00 

f-fia 

f;a0 

ffi 
J 

J a 
J 



f-af 

/,„ - (a - aHp)^J - HI 
t 1 

J,al3 - a 



'a/3 



{h - 2Ha + 2H(t> + </>) + {^H^i 



-f + 2af + af 

-f,a + a,af + Hf^a 



/ 



-/ + 2a/ + ^/ - (3if - d - ma + 3</) + ^x)/- 
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